function svl_jpe_figures_app()
% Appendix C

clear
clc
close all

[ell,theta,sigma,pi,rho,eta,i,M,A,B]=parameters();

% figure12(ell,theta,sigma,pi,eta,i,M,A,B)      % 280.645461 seconds
% figure13a(ell,theta,sigma,pi,eta,i,M,A,B)     %   1.287914 seconds
% figure13b(ell,theta,sigma,pi,rho,eta,i,M,A,B) %   1.305089 seconds

%=========================================================================
% parameters
%-------------------------------------------------------------------------
function [ell,theta,sigma,pi,rho,eta,i,M,A,B]=parameters()
% parameter setting
ell   = 0.5;       % probability of being C-type
theta = 0.5;       % C-type's bargaining power
sigma = 1;         % CRRA IES
pi    = 1;         % 1 - default rate
rho   = 0;         % recovery rate
eta   = 0;         % IRS
i     = 0.1; 
M     = 1; 
A     = 0.05;
B     = 0.05;

%=========================================================================
% subfunctions
%-------------------------------------------------------------------------
function u = u(q,sigma)
% CRRA utility: u(q)
if sigma < 1
    u = q.^(1-sigma)./(1-sigma);
elseif sigma == 1
    u = log(q);
end

function mu = mu(q,sigma)
% marginal utility: u'(q)
mu = q.^(-sigma);

function w = w(q,sigma,theta)
% w function in effective bargaining power
w = theta + (1-theta)*mu(q,sigma);

%=========================================================================
% solvers: figure 12
%-------------------------------------------------------------------------
function [G_eC,L_A,L_B,q0A,q1A,q0B,q1B,qBA,eN,dBA]=partialEquilibrium(eC,ell,theta,sigma,pi,eta,i,M,A,B)
% partial equilibrium, given eC

lb = [0,0,0,0,0,0,0]; % lower bound constraint
ub = [1,1,1,1,1,1,A/(1-eC)]; % upper bound constraint
x0 = [0.8,0.9,0.8,0.9,0.9,A/(A+B),A/(1-eC)/2]; % initial guess
opts = optimoptions(@fmincon,'Algorithm','sqp','Display','off'); 
if eC > 0 && eC < 1
    nonlcon  = @(x) fminconstr(x,eC,ell,theta,sigma,pi,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon,opts);
elseif eC == 0
    nonlcon0 = @(x) fminconstr0(x,eC,ell,theta,sigma,pi,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon0,opts);
elseif eC == 1
    nonlcon1 = @(x) fminconstr1(x,eC,ell,theta,sigma,pi,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon1,opts);
end

q0A = x(1);
q1A = x(2);
q0B = x(3);
q1B = x(4);
qBA = x(5);
eN  = x(6);
dBA = x(7);

if eC > 0 && eC < 1
    a_CA_n  = eN*(1-ell)/(eC*ell+eN*(1-ell))^(1-eta);
    a_CB_n  = (1-eN)*(1-ell)/((1-eC)*ell+(1-eN)*(1-ell))^(1-eta);
    a_CA_d  = 1-ell;
    S_CA    = theta*(u(q1A,sigma)-u(q0A,sigma)-q1A+q0A);
    S_CB_n  = theta*(u(q1B,sigma)-u(q0B,sigma)-q1B+q0B);
    S_CBA_d = theta*(u(qBA,sigma)-u(q0B,sigma)-qBA+q0B);
    L_A     = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
    L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1B,sigma,theta)*(mu(q1B,sigma)-1);
    S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
    S_CB_tilde  = - i*q0B - L_A*((1-theta)*(u(qBA,sigma)-u(q0B,sigma))+theta*(qBA-q0B)) - L_B_rho_bar*((1-theta)*(u(q1B,sigma)-u(q0B,sigma))+theta*(q1B-q0B)) ...
                  + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CA_d*S_CBA_d);
    G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / pi;
elseif eC == 0
    a_CA_n  = 0;
    a_CB_n  = 1-ell;
    a_CA_d  = 1-ell;
    S_CA    = theta*(u(q1A,sigma)-u(q0A,sigma)-q1A+q0A);
    S_CB_n  = theta*(u(q1B,sigma)-u(q0B,sigma)-q1B+q0B);
    S_CBA_d = theta*(u(qBA,sigma)-u(q0B,sigma)-qBA+q0B);
    L_A     = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
    L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1B,sigma,theta)*(mu(q1B,sigma)-1);
    S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
    S_CB_tilde  = - i*q0B - L_A*((1-theta)*(u(qBA,sigma)-u(q0B,sigma))+theta*(qBA-q0B)) - L_B_rho_bar*((1-theta)*(u(q1B,sigma)-u(q0B,sigma))+theta*(q1B-q0B)) ...
                  + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CA_d*S_CBA_d);
    G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / pi;
elseif eC == 1
    a_CA_n = 1-ell;
    a_CB_n = 0;
    a_CA_d = 1-ell;
    S_CA    = theta*(u(q1A,sigma)-u(q0A,sigma)-q1A+q0A);
    S_CB_n  = theta*(u(q1B,sigma)-u(q0B,sigma)-q1B+q0B);
    S_CBA_d = theta*(u(qBA,sigma)-u(q0B,sigma)-qBA+q0B);
    L_A     = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
    L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1B,sigma,theta)*(mu(q1B,sigma)-1);
    S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
    S_CB_tilde  = - i*q0B - L_A*((1-theta)*(u(qBA,sigma)-u(q0B,sigma))+theta*(qBA-q0B)) - L_B_rho_bar*((1-theta)*(u(q1B,sigma)-u(q0B,sigma))+theta*(q1B-q0B)) ...
                  + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CA_d*S_CBA_d);
    G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / pi;
end

function [eC,L_A,L_B,q0A,q1A,q0B,q1B,qBA,eN,dBA]=fixedPoint(ell,theta,sigma,pi,eta,i,M,A,B)
% general equilibrium 

% initial guess
eC = A/(A+B);

cnt = 0;
precision = 1;
while precision > 10^-6
    cnt = cnt + 1;
    if cnt > 1
        eC = temp_eC;
    end
    
    fprintf('finding a fixed point ... iter: %.0f ... ell = %.1f, theta = %.2f, pi = %.4f, eta = %.4f, i = %.4f, A= %.4f, B = %.4f ... eC = %.4f ',cnt,ell,theta,pi,eta,i,A,B,eC)
    [G_eC,L_A,L_B,q0A,q1A,q0B,q1B,qBA,eN,dBA]=partialEquilibrium(eC,ell,theta,sigma,pi,eta,i,M,A,B);
    fprintf('... G_eC = %.4f\n\n',G_eC)
    
    itr = 500; %1
    if cnt < itr
        if eC < 1
            if eC > 0
                temp_eC = max(0,min(1,eC + .5*G_eC)); 
%                 temp_eC = max(0,min(1,eC + 5*G_eC)); 
            elseif eC == 0
                break
            end
        elseif eC == 1
            break
        end
        precision = abs(G_eC);
    elseif cnt >= itr
        if cnt > itr
            temp_eC_2 = temp_eC_1;
        end
        if eC < 1
            if eC > 0
                temp_eC_1 = eC;
                temp_eC = max(0,min(1,eC + sign(G_eC)*10^-4)); 
            elseif eC == 0
                break
            end
        elseif eC == 1
            break
        end
        if cnt > itr
            if temp_eC_2 == temp_eC
                break
            end
        end
    end
end

%=========================================================================
% solvers: figure 13
%-------------------------------------------------------------------------
function S0=S0(sigma,ell,i)
% not participating
q  = (ell/(i+ell))^(1/sigma); % i = ell*(u'(q)-1)
S0 = - i*q + ell*(u(q,sigma)-q);

function [S1,S2]=partialEquilibrium_crs(e2,ell,theta,sigma,pi,eta,i,M,A,B)
% partial equilibrium, given e2

lb = [0,0,0,0, 0,0, 0,0,0, 0,0,0,0]; % lower bound constraint
ub = [1,1,1,1, 1,1, 1,1,1, Inf,Inf,Inf,Inf]; % upper bound constraint
x0 = [0.8,0.9,0.8,0.9, 0.5,0.5, 0.8,0.9,0.9, 2*A,2*B,A,B]; % initial guess
opts = optimoptions(@fmincon,'Algorithm','sqp','Display','off','FunctionTolerance',1e-10); %
nonlcon = @(x) fminconstr_crs(x,e2,ell,theta,sigma,pi,eta,i,M,A,B);
x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon,opts);

% segmented markets
q0A  = x(1);
q1A  = x(2);
q0B  = x(3);
q1Bn = x(4);
eC   = x(5);
eNn  = x(6);
eNd  = 1;
% consolidated market
q0   = x(7);
q1n  = x(8); 
q1d  = x(9); 
% d's
dA   = x(10);
dB   = x(11);
dA2  = x(12);
dB2  = x(13);
% phi
phi  = 1/M * ((1-e2)*(eC*q0A+(1-eC)*q0B)+e2*q0);

% matching probabilities
a_CA_n = (1-e2)^eta * eNn*(1-ell)/(eC*ell+eNn*(1-ell))^(1-eta);
a_CB_n = (1-e2)^eta * (1-eNn)*(1-ell)/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
a_NA_n = (1-e2)^eta * eC*ell/(eC*ell+eNn*(1-ell))^(1-eta);
a_NB_n = (1-e2)^eta * (1-eC)*ell/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
a_CA_d = (1-e2)^eta * eNd*(1-ell)/(eC*ell+eNd*(1-ell))^(1-eta);
a_NA_d = (1-e2)^eta * eC*ell/(eC*ell+eNd*(1-ell))^(1-eta);
a_C2   = e2^eta * (1-ell);
a_N2   = e2^eta * ell;

% segmented markets
L_A         = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1);
S_CA   = theta*(u(q1A, sigma)-u(q0A,sigma)-q1A +q0A);
S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
S_CA_tilde = - i*q0A - L_A*phi*dA + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
S_CB_tilde = - i*q0B - L_B_rho_bar*phi*dB + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n);
S_NA   = (1-theta)*(u(q1A, sigma)-u(q0A,sigma)-q1A +q0A);
S_NB_n = (1-theta)*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
S_NA_tilde = (1-ell)*pi*a_NA_n*S_NA   + (1-ell)*(1-pi)*a_NA_d*S_NA;
S_NB_tilde = (1-ell)*pi*a_NB_n*S_NB_n + (1-ell)*(1-pi)*a_NA_d*S_NA;
S1 = S_CA_tilde + S_NA_tilde;
% S1 = S_CB_tilde + S_NB_tilde;
S1 = S1 - S0(sigma,ell,i);

% consolidated market
L_A2         = ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1) + ell*(1-pi)*a_C2*theta/w(q1d,sigma,theta)*(mu(q1d,sigma)-1);
L_B2_rho_bar = ell*pi*a_C2*theta/w(q1n,sigma,theta)*(mu(q1n,sigma)-1);
S_C2_n = theta*(u(q1n,sigma)-u(q0,sigma)-q1n+q0);
S_C2_d = theta*(u(q1d,sigma)-u(q0,sigma)-q1d+q0);
S_C2_tilde = - i*q0 - L_A2*phi*dA2 - L_B2_rho_bar*phi*dB2 + ell*(u(q0,sigma)-q0+pi*a_C2*S_C2_n+(1-pi)*a_C2*S_C2_d);
S_N2_n = (1-theta)*(u(q1n,sigma)-u(q0,sigma)-q1n+q0);
S_N2_d = (1-theta)*(u(q1d,sigma)-u(q0,sigma)-q1d+q0);
S_N2_tilde = (1-ell)*a_N2*(pi*S_N2_n+(1-pi)*S_N2_d);
S2 = S_C2_tilde + S_N2_tilde;
S2 = S2 - S0(sigma,ell,i);

function [S1,S2]=partialEquilibrium_irs(e2,ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given e2

lb = [0,0,0,0,0,           0,0,0,       0,0,0,       0,0,0,0,         0,0,0,0,           ]; % lower bound constraint
ub = [1,1,1,1,1,           1,1,1,       1,1,1,       Inf,Inf,Inf,Inf, Inf,Inf,Inf,Inf,   ]; % upper bound constraint
x0 = [0.8,0.9,0.8,0.9,0.9, 0.5,0.5,0.5, 0.8,0.9,0.9, 2*A,2*B,A,B,     1e-2,1e-2,1e-2,1e-2]; % initial guess
opts = optimoptions(@fmincon,'Algorithm','sqp','Display','off','FunctionTolerance',1e-10);  %
nonlcon = @(x) fminconstr_irs(x,e2,ell,theta,sigma,pi,rho,eta,i,M,A,B);
x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon,opts);

% segmented markets
q0A  = x(1);
q1A  = x(2);
q0B  = x(3);
q1Bn = x(4);
q1Bd = x(5);
eC   = x(6);
eNn  = x(7);
eNd  = x(8);
% consolidated market
q0   = x(9);
q1n  = x(10); 
q1d  = x(11); 
% d's
dA   = x(12);
dB   = x(13);
dA2  = x(14);
dB2  = x(15);
% phi
phi  = 1/M * ((1-e2)*(eC*q0A+(1-eC)*q0B)+e2*q0);
% liquidity premia
L_A  = x(16);
L_A2 = x(17);
L_B_rho_bar  = x(18);
L_B2_rho_bar = x(19);

% matching probabilities
a_CA_n = (1-e2)^eta * eNn*(1-ell)/(eC*ell+eNn*(1-ell))^(1-eta);
a_CB_n = (1-e2)^eta * (1-eNn)*(1-ell)/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
a_NA_n = (1-e2)^eta * eC*ell/(eC*ell+eNn*(1-ell))^(1-eta);
a_NB_n = (1-e2)^eta * (1-eC)*ell/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
a_CA_d = (1-e2)^eta * eNd*(1-ell)/(eC*ell+eNd*(1-ell))^(1-eta);
a_CB_d = (1-e2)^eta * (1-eNd)*(1-ell)/((1-eC)*ell+(1-eNd)*(1-ell))^(1-eta);
a_NA_d = (1-e2)^eta * eC*ell/(eC*ell+eNd*(1-ell))^(1-eta);
a_NB_d = (1-e2)^eta * (1-eC)*ell/((1-eC)*ell+(1-eNd)*(1-ell))^(1-eta);
a_C2   = e2^eta * (1-ell);
a_N2   = e2^eta * ell;

% segmented markets
S_CA   = theta*(u(q1A, sigma)-u(q0A,sigma)-q1A +q0A);
S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
S_CB_d = theta*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);
S_CA_tilde = - i*q0A - L_A*phi*dA + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
S_CB_tilde = - i*q0B - L_B_rho_bar*phi*dB + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n+(1-pi)*a_CB_d*S_CB_d);
S_NA   = (1-theta)*(u(q1A, sigma)-u(q0A,sigma)-q1A +q0A);
S_NB_n = (1-theta)*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
S_NB_d = (1-theta)*(u(q1Bd,sigma)-u(q0B,sigma)-q1Bd+q0B);
S1 = S_CA_tilde + (1-ell)*pi*a_NA_n*S_NA + (1-ell)*(1-pi)*a_NA_d*S_NA;
% S1 = S_CB_tilde + (1-ell)*pi*a_NB_n*S_NB_n + (1-ell)*(1-pi)*a_NB_d*S_NB_d;
S1 = S1 - S0(sigma,ell,i);

% consolidated market
S_C2_n = theta*(u(q1n,sigma)-u(q0,sigma)-q1n+q0);
S_C2_d = theta*(u(q1d,sigma)-u(q0,sigma)-q1d+q0);
S_C2_tilde = - i*q0 - L_A2*phi*dA2 - L_B2_rho_bar*phi*dB2 + ell*(u(q0,sigma)-q0+pi*a_C2*S_C2_n+(1-pi)*a_C2*S_C2_d);
S_N2_n = (1-theta)*(u(q1n,sigma)-u(q0,sigma)-q1n+q0);
S_N2_d = (1-theta)*(u(q1d,sigma)-u(q0,sigma)-q1d+q0);
S2 = S_C2_tilde + (1-ell)*a_N2*(pi*S_N2_n+(1-pi)*S_N2_d);
S2 = S2 - S0(sigma,ell,i);

% liquidity premia
L_A  = 100 * L_A;
L_A2 = 100 * L_A2;
L_B  = 100 * L_B_rho_bar  / (pi+(1-pi)*rho);
L_B2 = 100 * L_B2_rho_bar / (pi+(1-pi)*rho);

%=========================================================================
% grid
%-------------------------------------------------------------------------
function [PI,EC,DBA]=PI_mat(ell,theta,sigma,pi,eta,i,M,A,B)

PI = [0 0.01 0.02 0.1:0.1:1];

EC  = zeros(1,length(PI));
DBA = zeros(1,length(PI));

for j=1:length(PI)
    pi = PI(j);
    
    [eC,L_A,L_B,q0A,q1A,q0B,q1B,qBA,eN,dBA]=fixedPoint(ell,theta,sigma,pi,eta,i,M,A,B);
    
    EC(j)  = eC;
    DBA(j) = dBA;
    
end

function [E2,S1_mat,S2_mat]=E2_mat_crs(ell,theta,sigma,pi,eta,i,M,A,B)
% partial equilibrium, given e2

E2 = 0:0.01:1;

S1_mat = zeros(1,length(E2));
S2_mat = zeros(1,length(E2));

for j=1:length(E2)
    e2 = E2(j);
    
    fprintf('Computing ... e2 = %.4f\n\n',e2)
    [S1,S2] = partialEquilibrium_crs(e2,ell,theta,sigma,pi,eta,i,M,A,B);
    
    S1_mat(j) = S1;
    S2_mat(j) = S2;
    
end

function [E2,S1_mat,S2_mat]=E2_mat_irs(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given e2

E2 = [0:0.047:1 1];

S1_mat = zeros(1,length(E2));
S2_mat = zeros(1,length(E2));

for j=1:length(E2)
    e2 = E2(j);
    
	fprintf('Computing ... e2 = %.4f\n\n',e2)
    [S1,S2] = partialEquilibrium_irs(e2,ell,theta,sigma,pi,rho,eta,i,M,A,B);

    S1_mat(j) = S1;
    S2_mat(j) = S2;

end

%=========================================================================
% figure 12
%-------------------------------------------------------------------------
function figure12(ell,theta,sigma,pi,eta,i,M,A,B)
figure()

A = 0.0109; B = 0.0477; 
pi = 0.996; rho = 0.4; sigma = 0.34; ell = 0.5; theta = 0.8; eta = 0.10; i = 0.0175;

[PI,EC,DBA]=PI_mat(ell,theta,sigma,pi,eta,i,M,A,B);
DA = (A-(1-EC).*DBA)./EC;

plot(PI,EC.*DA,'k--','Linewidth',3), hold on
plot(PI,(1-EC).*DBA,'k-','Linewidth',3), hold on
axis square
set(legend({'\,$$e_Cd_A$$','\,$$(1-e_C)d_{BA}$$'},'Interpreter','latex','Location','southeast','Orientation','vertical'),'FontSize',13)
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
xlabel('\fontsize{18} \pi'), 
set(gcf,'pos',[30 70 415 430])
% print('fig12','-depsc')

disp('done')

%=========================================================================
% figure 13
%-------------------------------------------------------------------------
function figure13a(ell,theta,sigma,pi,eta,i,M,A,B)
figure()

pi = 0.95; eta = 0;

[E2,S1_mat,S2_mat] = E2_mat_crs(ell,theta,sigma,pi,eta,i,M,A,B);

plot(E2,S1_mat,'k-','Linewidth',3), hold on
plot(E2,S2_mat,'k--','Linewidth',3), hold on
set(legend({'\,$$\mathcal{S}_1$$','\,$$\mathcal{S}_2$$'},'Interpreter','latex','Location','south','Orientation','vertical'),'FontSize',13)
yticks(linspace(1.885,1.925,9)*10^-3), yticklabels({' ',' ',' ',' ',' ',' ',' ',' ',' ',' '}), ylim([1.885,1.93]*10^-3)
axis square
set(gcf,'pos',[30 70 415 430])
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$e_2^{}$$','interpreter','latex'),'FontSize',15)
% print('fig13a','-depsc')

disp('done')

function figure13b(ell,theta,sigma,pi,rho,eta,i,M,A,B)
figure()

A = 0.0109; B = 0.0477;
pi = 0.996; rho = 0.4; sigma = 0.34; ell = 0.5; theta = 0.8; eta = 0.10; i = 0.0175; 

[E2,S1_mat,S2_mat]=E2_mat_irs(ell,theta,sigma,pi,rho,eta,i,M,A,B);

plot(E2,S1_mat,'k-','Linewidth',3), hold on
plot(E2,S2_mat,'k--','Linewidth',3), hold on
ylim([0 1.8]*10^-4), yticks([0:0.2:1.8]*10^-4), yticklabels({' ',' ',' ',' ',' ',' ',' ',' ',' ',' '})
set(legend({'\,$$\mathcal{S}_1$$','\,$$\mathcal{S}_2$$'},'Interpreter','latex','Location','south','Orientation','vertical'),'FontSize',13)
axis square
set(gcf,'pos',[480 70 415 430])
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$e_2^{}$$','interpreter','latex'),'FontSize',15)
% print('fig13b','-depsc')

disp('done')
